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Abstract 


In  this  report,  the  equations  of  motion  for  a  satellite 
in  restricted  four-body  motion  are  derived  and  examined.  An 
algorithm  is  developed  that  produces  a  minimum  AV  transfer 
trajectory  from  a  low  Earth  orbit  to  establishment  in  Wheeler's 
stable  periodic  orbit  around  L4.  Both  a  Hohmann  transfer  and 
an  infinite  velocity,  straight  line  transfer  are  examined  as 
initial  conditions  for  integration  of  the  equations  of  motion. 
Transfer  trajectories  are  determined  to  ten  points  around 
Wheeler's  orbit  so  that  a  minimum  AV  trajectory  can  be 
determined.  Two  types  of  minimum  AV  values  were  found. 

One  type  is  a  true  minimum  value,  the  other  type  is  a 
limiting  value  beyond  which  the  target  orbit  can  not  be 
reached.  Results  are  presented  on  a  transfer  trajectory 
that  gives  a  minimum  AV  from  the  points  examined.  From 
these  results  a  more  accurate  transfer  location  is  estimated. 
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DETERMINATION  OF  THE  MINIMUM  AV 
TRANSFER  TRAJECTORY  FROM  A  LOW 
EARTH  ORBIT  TO  A  STABLE  ORBIT 
AROUND  THE  LAGRANGIAN  POINT  L4 
IN  A  RESTRICTED  FOUR-BODY  SYSTEM 


I  Introduction 

Background 

In  the  past  few  years  the  idea  of  developing  space 
colonies  has  been  gaining  strength.  NASA  (Ref  6)  has 
recently  looked  at  the  feasibility  of  establishing  space 
settlements.  They  determined  that  space  colonization  is 
desirable  because  of  the  hope  it  offers  humanity.  It  has 
been  suggested  that  colonies  can  provide  a  solution  to  many 
earthly  problems  ranging  from  energy  sources  to  overcrowding 
population  growth.  The  location  of  such  colonies  has  been 
a  major  subarea  of  study.  The  majority  of  work  in  this  area 
has  been  to  find  periodic  orbits  in  the  restricted  three- 
body  problem,  with  the  Earth,  Moon  and  space  colony  forming 
the  three-body  model.  Recently,  however,  there  have  been 
attempts  to  investigate  the  motion  of  the  satellite  in  a 
restricted  four-body  model  by  adding  the  gravitational 
forces  of  the  Sun.  The  finding  of  periodic  orbits  in  the 
four-body  model  causes  questions  concerning  transfers  to 
these  orbits.  The  need  for  transporting  personnel,  mate¬ 
rials  and  equipment  to  the  site  at  an  efficient  and  economic 


rate  must  be  investigated.  That  is  the  main  purpose  of 
this  paper. 

From  a  military  viewpoint,  the  orbits  that  are  suitable 
for  nearby  colonies  may  also  be  likely  locations  for  military 
outposts  or  large  military  structures.  The  portion  of  the 
Earth  visible  at  any  one  time  from  the  orbit  and  its  prox¬ 
imity  to  botn  the  Earth  and  the  Moon  are  factors  that  would 
be  important  to  the  military  strategist.  For  these  reasons, 
orbits  in  the  vicinity  of  the  Lagrangian  points  L4  and  L5 
are  of  military  importance. 

Dr.  Gerard  K.  O'Neill  (Ref  8)  of  Princeton  University 
is  this  country's  leading  authority  in  the  field  of  space 
colonies  and  their  suitable  locations.  In  his  book  he 
states  that  L5  has  lost  its  status  as  the  probably  host  site 
for  the  first  space  colony.  However,  very  recently, 

J.  E.  Wheeler  (Ref  10)  investigated  the  existance  of  peri¬ 
odic,  stable  orbits  near  L4  and  L5.  His  finding  of  a 
stable  1/1  resonant  orbit  at  L4  revitalizes  the  concept  of 
a  colony  or  military  outpost  near  that  point.  This  1/1 
orbit  means  that  the  satellite  revolves  around  L4  once  for 
every  time  the  Moon  revolves  around  the  Earth. 

T.  A.  Heppenheimer  (Ref  5)  and  his  colleagues 
B.  O'Leary  and  D.  Kaplan  (Ref  7)  have  also  examined  the  pro¬ 
blem  of  colony  location  and  transfer  trajectories. 
Heppenheimer  found  a  2/1  resonant  orbit  around  the  Earth 
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that  could  be  reached  from  L2  by  a  Hohmann  transfer  sequence. 
He  did  not,  however,  examine  orbits  or  transfer  trajec¬ 
tories  in  the  neighborhood  of  L4  or  L5  because  current  mass- 
driver  capabilities  would  not  allow  a  catcher  near  either 
of  these  points.  In  their  work  together,  O'Leary,  Kaplan 
and  Heppenheimer  found  a  transfer  trajectory  from  L2  that 
circles  the  Moon  twice  and  then  passes  through  LI  to  enter 
a  2:1  resonant  orbit  around  the  Earth.  Again,  L4  and  L5 
were  not  examined.  As  the  accuracy  of  mass-driven  launchers 
improves,  orbits  for  manufacturing  facilities  near  L4  or  L5 
become  more  feasible.  As  this  happens,  the  need  for 
examining  transfer  trajectories  to  this  area  will  increase. 

As  the  need  for  examining  transfer  trajectories  grows, 
techniques  for  finding  these  orbits  must  be  developed.  Much 
of  today's  work  is  done  with  Hohmann  transfers,  such  as 
Heppenheimer ' s  work  (Ref  5).  Hohmann  transfers  provide 
the  minimum  AV  required  for  a  transfer  between  two  copla- 
nar,  circular  orbits  in  a  two-body  system.  In  a  restricted 
four-body  system,  such  as  exists  here,  Hohmann  transfers 
may  not  produce  the  minimum  AV  transfer.  The  necessity 
for  developing  a  better  method  definitely  exists. 
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Problem  and  Scope 

The  problem  is  to  determine  the  minimum  velocity  change 
required  to  transfer  from  a  low  Earth  orbit  to  the  orbit 
around  L5  found  by  Wheeler.  The  very  restricted  four-body 
model  was  used.  In  this  model,  the  Earth  and  Moon  travel 
in  circles  around  their  barycenter.  This  barycenter  and 
the  Sun  follow  circular  orbits  around  the  Earth-Moon-Sun 
barycenter.  The  motion  of  the  colony,  as  it  moves  through 
this  system,  is  governed  by  the  gravitational  forces  created 
by  each  of  the  other  three  bodies .  All  motion  in  the  system 
is  limited  to  one  plane. 

The  Lagrangian  Points,  L4  and  L5,  near  which  Wheeler 
found  his  orbit,  are  equilibrium  points  in  the  restricted 
three-body  problem.  At  these  points  the  sum  of  the  grav¬ 
itational  forces  of  the  Earth  and  Moon,  plus  any  centrif¬ 
ugal  force  on  the  satellite  equals  zero.  If  a  satellite  were 
placed  at  one  of  these  points  with  the  proper  initial  veloc¬ 
ity  it  will  remain  fixed  relative  to  the  Earth  and  Moon. 

When  the  gravitational  force  of  the  Sun  is  added  to  the 
model,  these  points  are  no  longer  equilibrium  points.  As 
shown  in  Fig  1,  the  points  L4  and  L5  form  a  equilateral 
triangle  with  the  Earth.  Therefore,  any  transfer  trajectory 
to  one  will  be  the  mirror  image  of  the  other.  The  only  time 
that  this  is  not  true  is  when  the  perturbing  effects  of 
other  heavenly  bodies  are  introduced.  These  perturbing 
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II  Problem  Analysis 


The  primary  background  work  that  this  thesis  is  based 
upon  was  done  by  J.  E.  Wheeler  (Ref  10)  as  an  AFIT  thesis. 
Wheeler's  dynamical  model  and  basic  problem  assumptions 
were  used  as  the  primary  foundations  for  this  paper.  From 
his  thesis  a  1/1  orbit  was  selected  as  a  target  for  a 
transfer  trajectory  from  a  low  Earth  orbit. 

Assumptions 

Following  are  the  basic  assumptions  used  in  the  very 
restricted  four-body  model. 

1)  The  Sun,  Earth,  and  Moon  are  considered  to  be  point 
masses. 

2)  The  mass  of  the  satellite  is  negligible  when  com¬ 
pared  to  the  other  three  bodies  and,  therefore,  it 
does  not  affect  their  motion. 

3)  The  gravitational  effects  of  the  other  planets  in 
our  solar  system  can  be  ignored. 

4)  The  motion  of  all  four  bodies  is  limited  to  one 
plane . 

5)  Earth  and  Moon  move  in  circular  orbits  about  their 
barycenter  at  a  constant  rotation  rate. 


6)  The  Earth-Moon  barycenter  moves  in  a  circular  orbit 
about  the  Earth-Moon-Sun  barycenter  at  a  constant 
rate. 

7)  The  satellite  moves  in  a  circular  orbit  about  the 
Earth  at  a  constant  rate. 

The  assumed  circular  orbits  are  actually  slightly 
elliptical.  The  Moon's  actual  eccentricity  is  0.0549 
(Ref  2:325),  and  the  Earth's  eccentricity  is  0.0167 
(Ref  2:360) .  The  slight  errors  that  result  from  assuming 
circular  orbits  will  be  on  the  order  of  the  eccentricities 
of  the  orbits.  The  assumption  of  planer  motion  is  also  not 
exact.  The  Moon's  orbit  is  inclined  5°8'  (Ref  2:325) 
from  the  Earth's  orbital  plane.  This  will  also  cause 
small  errors  to  occur.  These  errors  indicate  that  Wheeler's 
stable,  periodic  1/1  orbit  isn't  exactly  periodic  and 
minute  perturbations  will  exist.  These  small  perturbations 
will  also  be  present  in  any  transfer  trajectory  to  Wheeler's 
orbit  found  from  the  circular  restricted  four-body  model. 
Based  on  the  previously  mentioned  assumptions,  however, 
any  motion  found  from  the  circular  restricted  four  body 
model  should  be  nearly  identical  to  the  actual  motion. 
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Astronautical  Constants 


Before  continuing  with  the  development  of  the  equa¬ 
tions  of  motion,  the  parameters  used  in  the  analysis  roust 
be  defined.  In  Fig  2,  if  the  mass  of  the  Moon,  M 2,  is 
set  to  some  arbitrary  value,  y,  and  the  sum  of 
(Earth)  +  M2  is  equal  to  1,  then 


if  the  moment  arms  are  now  summed,  it  can  be  shown  that  the 
distance  from  the  Earth  and  the  Moon  to  the  barycenter  are 
y  and  1  -  y,  respectively.  The  solution  for  y  is 
then 

M 

11  "  M1  +^M2 

(2) 


The  value  of  y  used  is 
y  -  .0121396054 


(Ref  9:517) 


The  synodic  period  of  the  Moon,  T  ,  will  appear 

s 

throughout  this  report.  It  corresponds  to  9  which  is 
specified  as 


Pig  2.  Earth-Moon  System 
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with 


T  *  29.5305882  days 
s 


Now  referring  to  Fig  3, 


M1  f  M2  x 

M  328900.1; 

s 


M  ■  328900.12 
s 


(Ref  4:144) 


Te,  the  sidereal  year,  will  be  used  to  determine  &  , 


with 


T  *  365.25636556  days 
£< 


(Ref  3:334) 


since 


0  *  1  in  normalized  units,  then  a  becomes 


.  _  ,  2  tt  ,  ,29.5305882, 

a  "  '365.25636556'  *  1  2  ir  J 
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or 


d  =  .0808489351 

Again/  summing  moment  arms  in  Fig  3.  gives 
Rs  Ms  =  *B  (M1  +  V 

(6) 


or 


*B  MS 

KS  "  M1  +  M2 


328900.12 


Applying  Kepler ' s  Third  Law 

Tp  *  =  2  tt  (a3  /  GM)  ** 


or 


2  3 

<o  a  =  GM 


(7) 


(8) 


(9) 


The  rotation  rate  of  the  Earth-Moon  system  with  respect  to 
an  inertial  frame  is 


co  «  8  +  o 


(10) 
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Four-Body  Configuration 


So,  looking  at  just  the  Earth-Moon  system 


(6  +  ct)2  a  3  =  G  +  M2)  =  G  (1) 


With  a  equal  to  the  distance  from  the  Earth  to  the  Moon, 
1  Du,  then 


(1  +  d)2  =  G  =  1.16823442 


(ID 


Another  relationship  for  Rg  and  Rg  can  now  be 
solved  for  using  Kepler's  Law  again.  Concentrating 
(M1  +  M2)  at  Bem  and  rotating  about  BgEM  at  an  oi 
equal  to  &  gives 


(d2)  (Rg  +  Rg)3  —  G  (Mg  +  Mx  +  M2) 


(12) 


or 


(Rg  +  Rg)  =  RgE  =  388.8202849 


(13) 


and  from  Eq.  (7) , 


328900.12 
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Combining  Eq.s.  (7)  and  (13)  gives 


*  388.8191027 


Rs  =  .0011821799 


Equations  of  Motion 

All  motion  will  be  referred  to  the  inertial  reference 

frame  xT,yT  in  e  e  unit  vectors.  From  Fig  3, 

11  1'  2 

r  =  x  e  +  y  e 
1  2 


r  =  V  =  xe  +ye 
“  “  “1  2 


V  *  V  +  (0  x  r)  =  (x  -  8y)e  +  (y  +  0x)e 
“E  “  ~1  ~~2 


V  =  V  +  (a  x  R  )  =  [x  -  (9  +  a) y  +  a  sin  0]e 
~I  ~E  ~I  1 


+  [y  +  (9  +  a)x  +  a  Rg  cos  9]e 


The  kinetic  energy  is 


T  =  h  Mc  (Vr  .  Vx) 
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T  =  %  Mc  {[x  -  (0  +  a)y  +  a  Rg  sin  0 p 

•  •  •  •  M 

+  [y  +  (0  +  a)x  +  a  Rg  cos  0]^} 


The  potential  energy  is 

v  -  -  z  gm(Mc} 

i=l  ri 


V  =  -GM  { - --0—  P 

c  r  fv-n 1  2  +  v)2!3* 


[(x-y)  +y)  P  [  (1-y+x)  +  y  ] 


+  [(x+RgE  cos  0  +  (y-RgE  sin  0)  p} 


The  Lagrangian  is  L  =  T  -  V.  To  solve  for  the  equations 
of  motion,  Lagranges's  equations  are  employed: 


3L  d  ,  3L  .  n  •  i  •> 

‘  3t  (rr’  °-  1  1,2 


and  the  generalized  coordinates  are  q^  =  x,  q2  =  y 


d  L  d  3  jj 

5^  ~  (— )  =  0  yields  the  x  equation  of  motion 

9  X 


—  2y(8+a)  +  (0+a)2x  +  R_  a2  cos  8  -  — 


[  (x-y)2+y2]3/2 


Gy (1-y+x) 


GMS  (RSE  cos  8  +  x) 


[ (1-y+x) 2+y2] 3/2  [(RgE  cos  8+x) 2+(y-RSEsin8) 2]3/2 


Similarly, 


d_  ,3L. 
dt 


=  0  yields  the  y  equation  of 


motion 


-2  x  ( 8  +a )  +  ( 8  +a )  2y 


[  (1-y+x)2  +y2 


-  otRg  sin  8  - 


G(l-y)y 
[  ( x-y )  2+y2  ]3/2 


MsG  (y  -  Rse  sin  8) 

C  (RgE  cos  0+x)2  +  (y-RgE  sin  8)2  f^2 


To  put  these  equations  into  first-order  form  a  state 
space  approach  was  used.  The  state  vector  X  was  defined 
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The  following  first-order  differential  equations  results 


X  =  x  (25) 

1 

X  -  y  (26) 

2 

X  =  x  (27) 

3 

X  =  y  (28) 

4 


Note  that  from  Eq.  (3)  ,  d0  =  dt  or  6  =  t,  the  time,  in 
Eqs.  (22)  and  (23).  Therefore,  these  equations  vary  with 
time  and  may  be  represented  by  the  equation: 

x  (t)  =  f  [x  (t) ,  t] 

~  ~  (29) 

This  equation  says  that  the  first  order  non-linear  equations 
are  functions  of  the  states  themselves  and  of  time. 

For  small  displacements  from  the  states  at  any  time,  t, 
Eq.  (29)  becomes: 

x  (t)  +  6x  (t)  *  f  [(x(t)  +  6x  (t)  ),  t] 

(30) 

By  the  use  of  a  Taylor  series  expansion  and  Eq.  (29) , 

Eq.  (30)  becomes: 

si  (t>  “  ^fft)  (») 

(31) 
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This  equation  relates  small  displacements  in  the  states  to 
the  equations  of  motion  through  a  square  matrix  A(t),  where 


Alt)  =  Cr=  ]  -  [•!#] 


x  3*j 


(32) 


In  this  case  A  is  a  4x4  matrix  that  is  used  to 
linearize  small  displacements.  Its  complete  listing  can  be 
found  in  Appendix  A. 


Equations  of  Variation 

The  method  used  by  Wheeler  to  produce  periodic  motion 
is  almost  the  same  as  the  method  used  to  produce  a  transfer 
trajectory  to  a  specific  point.  In  order  to  accomplish 
this,  the  initial  conditions  must  be  related  to  the  state 
of  the  system  at  the  time,  t.  Starting  with  an  initial 
trajectory  xQ  (t) ,  the  states  can  be  perturbed  to  give  a 
general  solution  close  to  xQ  (t) . 

x  (t)  =  Xq  (t)  +  6  x  (t) 

(33) 


Expanding  in  a  Taylor  series  gives 

x  (t)  -  x  (t)  +  I  xft  <5x(t,.)  +  higher  order  terms 

~°  9£(t0)l  0  0 


Utilizing  the  state  transition  matrix,  $  (t,  tp)  which  is 
defined  by 


/ 


n 


6x(t0)  =  t  (t,  tQ)  6x(tQ) 


(34) 


or 


*(t'V 


SXjt) 

3£<V 


(35) 


This  matrix  relates  changes  in  the  states  at  t  to  changes 
in  the  initial  states  at  t^.  Differentiating  Eq.  (34) 
yields. 

6x(t)  =  $(t,tQ) 6xQ (tQ)  +  $(t,tQ)  6x (tp) 
or 

6x(t)  =  $(t,t  )6x  (t  ) 

~  0  0 

(36) 

since  <5x(t0)  equals  zero.  Placing  Eqs.  (34)  and  (36)  into 
Eq.  (31)  gives 

4(t,tQ)  =  A(t)  $(t,tQ) 

(37) 

The  last  equation  allows  the  propogation  of  4>  along  the 
trajectery  to  be  determined. 

Simultaneous  numerical  integration  of  the  equations  of 
motion  and  variation  will  provide  the  needed  information  on 
a  particular  solution  at  the  final  time,  so  that  a  transfer 
trajectory  can  be  produced. 


19 


Ill  Transfer  Trajector 


The  problem  of  finding  an  optimum  transfer  trajectory 
to  Wheeler's  orbit  can  be  divided  into  three  sub-problems. 
The  first  part  is  a  two-point  boundary  value  problem 
attempting  to  arrive  at  a  designated  point  on  Wheeler's 
orbit  at  a  specified  time.  The  second  part  of  the  problem 
consists  of  minimizing  the  delta  V  for  the  transfer  traj¬ 
ectory  found  in  the  two-point  boundary  value  problem.  The 
final  part  of  the  problem  is  to  repeat  parts  one  and  two  at 
representative  points  around  Wheeler's  orbit.  From  the 
collection  of  transfer  trajectories  the  minimum  cost  trans¬ 
fer  to  Wheeler's  orbit  can  be  determined. 

Two-Point  Boundary  Value  Problem 

This  two-point  boundary  value  problem  is  bounded  by  a 
known  starting  point  r(tQ)  and  a  known  ending  point 
r(t1).  The  orbit  velocity,  v(t1),  is  also  known.  The 
velocity  at  launch,  v(t^),  needs  to  be  found  such  that  a 
transfer  is  made  from  r(tQ)  to  rft^.  There  are  three 
launch  parameters  that  can  be  specified  to  help  determine 
the  initial  conditions  for  launch.  See  Fig  4.  The  first 
parameter  is  the  impulsive  velocity  change,  that  is 

added  to  the  satellites  orbital  velocity  around  the  Earth 
to  cause  it  to  leave  its  low  Earth  orbit.  The  second 
parameter  is  the  launch  angle,  that  determines  where 
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juju  y  wm  jpu 


in  the  orbit  the  impulse  should  be  added.  The  final  param¬ 
eter  is  the  launch  time,  t^.  This  time  tells  where  the 
Sun  is  in  relation  to  the  other  three  bodies.  Wheeler's 
orbit  is  phase  dependent,  so  the  position  of  the  influencing 
bodies  is  very  important.  The  satellite  must  arrive  at  a 
point  on  Wheeler's  orbit  when  the  Earth,  Moon  and  Sun  are 
in  the  required  positions  for  that  point  on  the  orbit.  That 
makes  the  time  of  arrival,  tQ,  very  important,  but  for  a 
specified  AV  and  <f>  a  specific  tQ  is  required  to  arrive 
at  the  target  at  t^ . 


Fig  4.  Launch  parameters 
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The  method  used  to  produce  a  transfer  trajectory  that 
solves  the  two-point  boundary  value  problem  is  very  similar 
to  that  used  by  Wheeler  to  find  his  periodic  orbit.  The 
end  point  is  not  the  same  as  the  starting  point,  as  in 
Wheeler's  problem,  but  the  concept  is  the  same.  By  start¬ 
ing  with  an  estimate  of  the  initial  state,  x (t^) ,  the 
states  can  be  propogated  along  the  trajectory  by  Eq.  (29) . 

x(t)  =  f  [x(t) ,  t] 

At  x(t^),  it  is  found  that  the  trajectory  does  not  close 
on  the  desired  point,  x(t^),  but  is  some  Ax  distance 
away  from  the  target. 

The  error  vector  is  then  defined  by 

Ax  =  5(^1  -  xft^ 

(38) 

To  close  the  trajectory  on  x(t^).  Ax  must  be  driven  to 
zero.  This  can  be  done  by  adding  a  6x(tg)  which  will 
generate  a  6x(t1)  as  expressed  by  Eq.  (34) 

fix(t1)  =  <Sx(to) 

So,  to  obtain  closure. 


Ax  *  5x(t^) 


(39) 
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or 


Ax  =  $(t1,tQ)  5x(tQ) 

(40) 

However,  the  initial  state,  x(tQ)  ,  is  dependent  upon  the 
launch  parameters,  L.  Therefore, 

5x(tQ)  =  C5L 

(41) 

Where  C  is  a  2x4  matrix  relating  changes  in  the 
launch  parameters  and  changes  in  the  initial  state. 

So, 

Ax  *  $C6L 

(42) 


or 


5L  =  B-1Ax 


(43) 


where 

B  =  <&C. 

(44) 

From  this  last  equation,  corrections  to  the  launch 
parameters  can  be  solved  for  in  terms  of  known  quantities. 
An  algorithm  for  solving  this  two-point  boundary  value 
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problem  and  producing  a  transfer  trajectory  is  summarized 
as  follows: 

Equations  of  motion:  x( tQ)  =  f  [x(t  )  ,t] 

where  x(t^)  equals  initial  conditions. 

Equations  of  variation:  $  =  A4> ,  $  (tQ  ,tQ)  =  I 

At  t  =  0,  a  guess  is  made  of  the  launch  parameters,  L, 
which  determine  the  initial  conditions,  x(t^).  Integrating 
the  equations  of  motion  and  the  equations  of  variation  to 
some  time,  t  =  t^,  the  resulting  quantities  are  xft^) 
and  <J(t^,tQ).  Utilizing  Eq.  (44),  B  can  be  found.  Also, 
knowing  that  Ax  =  xft^)  -  x(t^) ,  6L  can  be  computed  from 
Eq.  (43) .  The  next  estimate  of  the  launch  parameters  is 
then 

*=  L  -  SL 

(45) 

which  will  give  us  x^tg).  This  algorithm  will  then  repeat 
itself  until  Ax  is  smaller  than  some  convergence  criteria, 
thus  defining  a  transfer  trajectory.  Appendix  B  has  a 
complete  development  of  the  B  matrix  that  was  used. 
Appendix  C  is  a  computer  listing  of  this  algorithm. 

Minimum  AV 

The  second  part  of  the  problem  is  to  minimize  the  total 
velocity  change  required  to  establish  the  satellite  in 
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•  Wheeler's  orbit  at  a  particular  target  point.  The  object 


of  this  sub-problem  is  displayed  in  Pig  5. 


Fig  5.  Minimize  total  AV 

This  figure  displays  the  type  of  results  that  are  ex¬ 
pected.  As  the  TOF  approaches  zero  on  the  graph,  the 
required  velocity  to  get  to  the  point  would  have  to  increase. 
The  limiting  case  would  be  for  an  infinite  velocity,  which 
would  result  from  an  infinite  AV.  An  infinite  velocity 
would  produce  a  straight-line  transfer  and  it  would  happen 
instantaneously.  So,  as  TOF  approaches  zero,  AV  would 
approach  infinity. 

The  minimum  AV  point  on  the  graph  may  be  character¬ 
ized  by  two  types  of  behavior.  It  may  truly  be  a  minimum 


value  and  past  the  corresponding  TOF  AV  would  increase. 
(Curve  A).  However,  there  is  another  possibility.  The 
minimum  AV  may  be  a  limiting  value.  Any  TOF  beyond 
this  point  would  produce  a  lower  AV,  but  the  resulting 
velocity  would  not  be  sufficient  to  reach  the  target. 

(Curve  B) .  A  Hohmann  transfer  is  an  example  of  this  in  two 
body  motion.  It  provides  the  minimum  AV  trajectory  that 
will  complete  the  transfer  (Fig  6) .  Anything  less  falls 
short  of  its  target  (Fig  7) . 

To  produce  a  curve  similar  to  Fig  5,  small  changes  in 
the  TOF  used  in  the  algorithm  that  solves  the  two-point 
boundary  value  problem  will  be  made  to  produce  changes  in 
the  total  AV.  This  process  is  repeated  until  the  minimum 
is  bracketed  or  it  becomes  apparent  that  a  limiting  value 
exists. 

If  TOF ( 1 ) <  TOF ( 2 ) <  TOF ( 3 ) ,  AV(1)>  AV(2)  and  AV(2)< 

AV(3),  then  a  minimum  exists  between  T0F(1)  and  T0F(3). 
From  this  point  a  bisection  routine  is  used  to  find  the 
minimum  AV.  The  bisection  algorithm  is  as  follows: 

1)  TOF  =  ( TOF { 1 )  +  TOF ( 3) ) /2 

2)  Use  algorithm  for  boundary  value  problem  to  find 
AV  for  TOF. 

3)  This  TOF  replaces  T0F(1)  or  T0F(3),  whichever  has 
the  largest  associated  AV.  e.g.  If  AV(1)<  AV(3) ,  then 
TOF  -*■  TOF (3)  and  the  algorithm  is  repeated. 
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4)  This  procedure  will  continue  until  the  difference 
between  the  existing  minimum  AV  and  the  AV  for  TOF 
is  less  than  a  specified  convergence  criteria. 

If  a  limiting  value  appears  to  exist,  the  original 
procedure  of  increasing  TOF  by  small  amounts  is  continued 
until  the  target  point  can  no  longer  be  reached  to  the 
desired  accuracy. 

Minimum  transfer 

After  finding  the  minimum  transfer  trajectory  to  one 
point  on  Wheeler's  orbit,  the  process  is  repeated  for  a 
nearby  point.  This  method  continues  until  the  entire  orbit 
has  been  examined.  It  must  be  remembered  that  Wheeler's 
orbit  is  phase  dependent,  and  therefore  so  are  the  transfer 
trajectories.  That  is  why  the  entire  orbit  must  be  examined. 
By  utilizing  only  a  representative  number  of  points  around 
the  orbit  a  curve  relating  AV  to  time  of  arrival  (position) 
on  the  orbit  can  be  developed.  See  Fig  8.  For  this  curve, 
minimum  points  can  be  examined  further  for  more  precise 
results.  Because  Wheeler's  orbit  is  periodic,  the  points 
at  each  end  of  the  curve  will  be  the  same. 
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Fig  8.  Minimum  &V 
Initial  Conditions 

To  begin  integration  for  this  problem,  approximate 
initial  conditions  must  be  determined.  Launch  conditions 
must  be  established  that  will  allow  the  satellite  to  arrive 
at  the  target  point  when  the  Earth,  Moon  and  Sun  are  in  the 
proper  position  to  allow  maneuvering  into  the  orbit.  This 
can  be  done  by  specifying  the  values  of  0  and  a.  See 
Fig  3.  However,  0  =  t  (Ref  10:19)  and,  by  the  same 
reasoning,  specifying  time  will  also  determine  a.  It  is, 
therefore,  necessary  that  the  values  of  0  and  a  at  t  *=  0 
be  known  to  correctly  relate  the  four  bodies.  Initial 
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conditions  on  the  satellite  must  also  be  determined.  Its 
orbital  altitude  and  velocity  are  necessities  for  formu¬ 
lation  of  the  proper  launch  conditions.  The  final  initial 
conditions  are  the  launch  parameters  themselves.  The 
launch  angle,  <p ,  and  impulsive  velocity  change,  AV^ 
are  required  for  integration  of  the  transfer  trajectory. 

Because  Wheeler's  orbit  is  phase  dependent,  the  posi¬ 
tions  of  the  Earth,  Moon  and  Sun  determine  when  a  point  is 
on  the  orbit  and  when  it  is  not.  Therefore,  it  is  imper¬ 
ative  that  the  satellite  arrive  at  the  target  point  when 
the  three  governing  bodies,  (Earth,  Moon  and  Sun),  are  in 
positions  that  allow  that  point  to  be  part  of  Wheeler's 
periodic  stable  orbit.  The  time  of  arrival,  t^,  is  re¬ 
lated  to  the  launch  time,  t^,  by  the  time  of  flight,  TOF. 

TOF  =  t  -  t 
1  0 

(46) 

To  ensure  compatability  of  time  measurements,  Wheeler's 
initial  conditions  of  0  =  a  =  0  at  t  =  0  were  used. 

They  are  diagrammed  in  Fig  9.  From  this  base,  t^,  which 
expresses  the  Earth,  Moon,  Sun  positioning,  can  be  deter¬ 
mined  for  any  point  on  the  target  orbit.  Placing  a  spec¬ 
ified  TOF  and  t^  into  Eq.  (46)  yields  a  launch  time, 
tQ.  This  value  can  then  be  used  in  the  equations  of  motion 
(Eqs.  25-28) ,  for  the  integration.  The  target  orbit  and 


Fig  9.  Initial  Condition  Configuration 
of  Four-Body  Model 
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t^'s  for  selected  points  are  shown  in  Pig  10  and  Table  1. 
These  were  generated  from  an  adaptation  of  Wheeler's  pro¬ 
gram  which  is  listed  in  Appendix  C. 

Four  conditions  must  be  established  for  the  satellite 
before  its  launch  position  can  be  determined.  Looking  at 
Eqs.  (B-l,  2,  7,  8),  it  is  apparent  that  the  orbital 
altitude  “*1,  orbital  velocity  -Vc,  launch  angle  -<J>, 
and  velocity  change  -AV,  are  necessary  to  compute  the 
initial  conditions  of  the  satellite  for  the  transfer  inte¬ 
gration.  The  orbital  altitude  selected  was  r^  =  .0170110431 
Du  in  Wheeler's  units  where  1  Du  is  the  Earth-Moon  dis¬ 
tance  of  384,400  Km  (1:441).  This  corresponds  to  a  100 
mile  orbit  which  is  a  commonly  used  low  Earth  orbit. 

The  orbital  velocity  of  the  satellite  is  dependent  upon 
the  altitude  chosen.  Assuming  two-body  motion  because  of 
the  satellites  proximity  to  the  Earth,  the  orbital  velocity 
is  given  by 

Vc  -  (Ref  2:165) 

(47) 

where  the  Earth's  gravitational  parameter  is  defined  as 

Ue  =  GMX  (Ref  2:14) 

(48) 
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Target  Points  On  Orbit  Around  L4 


/ 


with  G  =  1.16823442 

(ID 

and  Mx  *  1  -  y  =  .985818095 

(1) 

then,  =  1.15166663 

So,  placing  these  values  in  Eq.  (47) 

Vc  =  8.22806889  Du/Tu 

Even  if  a  four  body  model  is  used,  the  close  Earth  region  is 
dominated  by  the  Earth,  and  the  optimum  transfer  orbit  should 
not  depend  strongly  on  Vc. 

With  Vc  and  r^,  determined,  it  is  now  necessary  to  come 
up  with  a  good  initial  guess  for  <j>,  AV^  and  TOF.  Two 
possible  methods  may  be  used.  One  method  is  to  ignore  the 
effects  of  the  Moon  and  the  Sun,  making  the  transfer  a  two 
body  problem.  Under  these  conditions  a  Hohmann  transfer 
could  be  used  tc  provide  the  initial  <t>,  AV^,  and  also  the 
TOF.  A  minimum  AV^  Hohmann  transfer  in  a  non- rotating 
frame  would  require  that  the  satellite  be  launched  from  low 
Earth  orbit  when  it  preceded  the  target  by  180°.  See  Fig  10. 
Using  point  80  from  Fig  10  and  Table  1, 

x  =  -.6062887269693 
y  -  .6535319881417 
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Prom  Fig  11/ 


tan  a  =  y/(x  +  y) 

a  =  46.58° 

<j>  =  360  -  a  =  313.42° 


(48) 


(49) 


or 


<j>  =  5.470196073  radians 


However,  this  is  a  rotating  reference  frame.  Once  the  TOF 
is  determined  a  correction  can  be  added  to  <f>  for  the  amount 
of  rotation  during  the  travel  of  the  satellite. 


TOF  is  determined  by  utilizing 
TOF  =  n  (a ^  /  u)* 


(Ref  2:165) 

(50) 


at  =  (r  +  r  )/2  =  .458382633 
^  1  2 

TOF  =  .908507911  Tu 


(Ref  2:164) 

(51) 


This  is  approximately  4.26993184  days,  as  1  Tu  equals 
4.69993908  days  (Ref  10:75).  This  is  based  upon  the  fact 
that  the  period  of  the  stable  orbit  is  the  same  as  the 
synodic  period  of  the  Moon  at  2tt  Tu's. 


/ 


Using  the  TOF,  a  correction  to  <f>  may  now  be  calcu¬ 
lated.  During  .908507911  Tu's  the  system  will  rotate 
.908507911  radians.  This  value  should  be  added  to  the  <J> 
for  no  rotation  to  give  a  corrected  value, 

<j>  =  5.470196073  +  .908507911 

<t>  =  .095518677  radians 

Because  the  total  energy  of  the  Hohmann  transfer  is 
constant, 

Et  =  /  (rx  +  r2)  (Ref  2:164) 

(52) 

rx  *  1.70110431  x  10-2  Du 

r2  ■  (  (x  +  y)2  +  y2  )*  m  .899754223  Du 

Ue  =  1.15166663 
Then, 

Et  =  -1.256228473 

The  velocity  at  the  launch  point  of  the  transfer  tra¬ 
jectory  can  now  be  found. 

vl=  C2(y@/r  +  Et)]J*  (Ref  2:164) 

(53) 
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r 


V  =  11.52778289  Du/Tu 


Then, 


V  =  3.299713999  Du/Tu 
c 


The  other  technique  used  to  find  approximate  initial 
conditions  was  to  assume  a  straight-line  transfer  orbit. 

See  Fig  12.  This  type  of  trajectory  means  that  <j>  would  be 
determined  by  where  the  straight-line  transfer  was  tangent 
to  the  low  Earth  orbit.  From  Fig  12  it  can  be  seen  that 


sin  0  =  r  /r 
1  2 

8  =  1.083317004'’ 


(54) 


tan  a  =  (y/x  +  \i) 

a  =  46.58085192 
<t>  -  90  -  a  +  0 

<p  =  44.5024°  =  .776714541  rad 


(55) 


(56) 


To  actually  achieve  a  straight-line  transfer  trajectory 
on  infinite  velocity  would  be  necessary.  To  represent  a 
infinite  velocity,  Vfc  =  1  x  105  Du/Tu' s  was  chosen.  There¬ 
fore,  “lx  105  Du/Tu.  To  find  TOF, 
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TOF  =  r  /V 
3  t 


-  U22  -  r!2)^Vt 


=  8.99690466  x  10"6  Tu 


(57) 


Methodolgy 

The  initial  approach  to  the  problem  began  by  using  a 
Hohmann  transfer  to  approximate  the  initial  transfer  tra¬ 
jectory.  The  first  integration  was  run  with  an  accuracy  of 

10-1  Du.  From  these  results,  another  integration  was  run 

—2 

with  an  accuracy  of  10  Du.  However,  this  accuracy  could 
not  be  achieved.  Plotting  the  transfer  trajectory  showed 
that  it  never  approached  target  point  80  but  appeared  to  be 
drawn  away  by  the  Sun's  gravitational  attraction.  This 
would  indicate  that  the  Sun's  perturbations  are  too  large 
to  ignore,  even  for  determining  an  initial  guess  at  a  trans¬ 
fer  trajectory. 

A  Hohmann  transfer  is  not  a  member  of  the  solution  set 
to  the  four-body  equations  of  motion  because  of  the  assump¬ 
tion  that  the  Sun's  effect  can  be  ignored.  However,  a 
straight-line  trajectory  is  an  actual  solution,  even  if  it 
can  not  be  practically  achieved.  Therefore,  a  straight-line 
trajectory  was  assumed  and  the  algorithm  was  run  for  target 
point  80. 
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The  first  run  for  a  straight-line  transfer  was  run  to 
10  ®  Du  accuracy  which  is  less  than  a  quarter  mile.  This 
integration  verified  that  the  algorithm  was  at  least  func¬ 
tioning  properly  as  it  produced  a  transfer  trajectory. 

Having  found  a  transfer  to  point  80,  it  was  necessary  to 
step  down  AV  until  a  minimum  could  be  located.  AV  was 
stepped  down  from  10^  to  10^,  then  to  10^  with  acceptable 
results.  For  each  AV  change  a  new  integration  time  was 
computed.  Throughout  these  integrations  the  initial  <J> 
value  used  remained  the  same.  At  AV  =  10  the  accuracy 
was  increased  to  10  and  the  procedure  of  using  the 
results  of  one  integration  as  input  for  the  next  was  begun. 
It  was  also  at  this  point  that  the  integration  time  or  TOF 
became  the  primary  control  variable.  By  extrapolation  from 
previous  results,  a  new  set  of  AV.^  and  <J>  were  formulated 
for  a  specified  TOF.  This  procedure  was  continued  until 
the  minimum  AV  was  bracketed.  This  was  indicated  by  an 
increase  in  AV  from  an  increase  in  TOF.  From  this  point 
on,  a  bisection  method  was  utilized. 

It  was  during  the  use  of  extrapolation  procedures  that 
accuracy  became  a  problem.  When  using  the  extrapolated 
values  as  initial  integration  conditions  it  was  sometimes 
necessary  to  drop  the  desired  accuracy  to  10"1  Du  just  to 
get  some  kind  of  a  transfer  orbit  for  results.  Using  the 
resulting  and  AV  from  this  orbit  the  accuracy  could 
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be  increased  by  one  or  two  orders  of  magnitude  until  10“^ 
was  achieved.  At  10  5  Du  accuracy  (about  2.5  miles)  a 
new,  larger  TOF  would  be  chosen  and  a  new  <J>  and  AV^ 
values  computed. 

The  bisection  routine  used  consisted  of  choosing  a  TOF 
midway  between  the  two  bracketing  values.  When  the  result¬ 
ing  AV  was  formulated  it  replaced  the  highest  of  the  two 
bracketing  AV's.  This  process  was  then  repeated  until  the 

difference  between  the  bracketing  AV's  was  less  than 
—  9 

10  Du/Tu. 

During  the  process  of  finding  transfer  trajectories  to 
points  10,  20,  30  and  100  it  became  apparent  that  the 
AV  vs  TOF  curve  did  not  have  a  minimum,  but  actually 
appeared  to  approach  a  limiting  value  as  was  discussed 
previously.  The  primary  indication  was  the  inability  to 
maintain  a  nonsingular  B  matrix  past  a  certain  TOF  for 
the  desired  level  of  accuracy.  What  this  means  is  that  past 
a  certain  TOF  the  AV's  that  are  found  are  not  large 
enough  to  propel  the  satellite  to  within  the  required  accu¬ 
racy  distance  of  the  target  point.  Since  a  bisection  routine 
can  not  be  used  here,  a  new  technique  was  required.  When  it 
became  obvious  that  a  minimum  was  not  present  and  a  limiting 
value  was  present,  the  last  TOF  that  produced  an  accuracy 
within  10  Du  was  noted.  To  this  TOF  decreasing  values 
of  ATOF  were  added  until  10  ®  accuracy  could  be  attained. 
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This  procedure  was  continued  until  a  ATOF  of  10 
duced  a  AV  change  of  less  than  10  .  As  an  example: 

1.  Last  TOF  with  10“9  accuracy  -  1.6  Tu 

2.  Add  .01  Tu  to  get  next  attempted  value  -  1.61  Tu 

3.  Extrapolate  from  previous  findings  to  get  <j>  and  AV^ 

_  Q 

4.  Attempt  to  integrate  to  10  accuracy 

5.  If  successful,  repeat  1  -*■  4 

6.  If  not  successful,  change  ATOF  to  .001  -  1.601  Tu 

7.  Repeat  3  and  4 

8.  If  successful  repeat  at  that  level  until  10  can 
not  be  reached,  (i.e.  -  1.607  good,  1.608  not  good) 

9.  Then  proceed  to  1.6071 

-9 

10.  Continue  process  until  ATOF  is  10  and  the  change 
in  AV  <  10~9 

After  the  minimum  AV  to  each  arrival  point  has  been 
found,  a  graph  can  be  constructed  that  will  allow  an  esti¬ 
mate  of  the  arrival  point  with  the  minimum  AV.  (See  Fig  8) 
If  the  apparent  minimum  is  not  at  one  of  the  examined  points, 
several  points  near  the  estimated  minimum  may  be  examined  to 
arrive  at  the  transfer  trajectory  with  a  minimum  AV.  It 
must  be  remembered  that  between  any  two  points  on  the  target 
orbit  there  can  be  found  another  point,  so  this  procedure 
could  continue  forever.  However,  the  target  orbit  was  des¬ 
cribed  by  100  points  and  10  points  were  examined.  That  will 
be  the  limiting  spacing  for  this  thesis. 
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IV  Results  and  Discussion 


Results 

A  transfer  trajectory  was  found  to  each  of  the  ten 
points  examined  on  Wheeler's  orbit.  Each  of  these  points 
represents  an  arrival  time  and  Table  1  shows  these  rela¬ 
tionships.  All  ten  transfer  trajectories  can  be  character¬ 
ized  by  a  unique  set  of  launch  parameters  and  maneuver 
characteristics.  These  are  tabulated  in  Table  2  for  easy 
reference.  Of  the  ten  orbits,  the  transfer  to  point  10  has 
the  lowest  total  Av. 

AV  =  3.972123159  Du/Tu  =  3760.10622  m/sec 
Its  launch  parameters  are 

<J>  =  .724843243  radians  or  41.5204586° 

AVX  =  3.305668879  Du/Tu  or  3129.22475  m/sec 
The  time  from  launch  to  arrival  at  Wheeler's  orbit  is 

TOF  =  1.80201  Tu  or  8.46933723  days 
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To  enter  Wheeler's  orbit  after  arriving  at  point  10  a 
second  burn  is  required.  It  will  produce  AV ,  which 
is 


AV2  =  .666454280  Du/Tu  =  630.881  m/sec 

All  ten  transfer  trajectories  are  shown  in  Figs  (13-22) . 

In  these  figures  the  Earth  is  located  at  (0.121396054,  0.0) 
and  the  Moon  at  (-0.9878603946,  0.0). 

Discussion 

If  the  information  found  in  Table  2  is  plotted  on 
curves,  it  is  easier  to  establish  relationships  between 
orbits,  or  at  least  to  attempt  to  understand  why  some  of 
the  results  are  as  they  are. 
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Minimum  OeltaV  Transfer  To  Point  40 (T»2. 513274122) 

Pig  16. 
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Minimum  DeltaV  Transfer  To  Point  60 (T=3. 76991118) 

Fig  18. 
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Minimum  DeltaV  Transfer  To  Point  90 (T=5. 654866776) 


Fig  21 
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TRANSFER  TRAJECTORY 
TO  ORBIT 


Looking  first  at  the  plot  of  total  AV  versus  the 

target  point  (Fig  23)  and  comparing  it  to  plots  of  its 

components,  AV^,  and  AV2  (Figs  24,  25),  it  is  possible 

to  see  which  of  the  two  burns  was  most  influencial  in 

determining  the  minimum  total  AV.  As  can  be  seen,  the 

total  AV  curve  is  in  nearly  perfect  phase  agreement  with 

the  AV  curve.  AV  ,  however,  is  almost  180°  out  of 
2  1 

phase.  This  would  indicate  that  AV^  was  predominant  in 
determining  the  shape  of  the  total  AV  curve,  even  though 
the  AV^  only  ranges  from  20-27  percent  of  AV^  in  all 
cases.  This  can  be  explained  by  looking  at  Wheeler's 
periodic  orbit  and  examining  the  target  points  (Fig  10) . 
Those  points  requiring  the  largest  AV^,  points  30  and  40, 
are  the  farthest  from  the  Earth,  while  those  with  the  lowest 
AV^,  points  60  and  70,  are  the  closest.  However,  because 
of  the  clockwise  motion  of  Wheeler's  orbit,  the  velocity 
required  to  enter  the  stable  orbit  on  the  Earth  side  is 
nearly  opposite  that  of  the  arriving  satellite  and  a  large 
AV2  is  necessary.  In  contrast,  the  velocities  at  the  more 
distant  points  on  the  far  side  of  the  orbit  are  more  closely 
aligned  with  the  incoming  satellite's  velocity  and  a 
smaller  AV2  results.  This  explaines  why  AV^  and  AV2 
are  out  of  phase  so  dramatically  and  demonstrates  the  pre¬ 
dominance  of  AV2  in  shaping  the  total  AV  curve. 
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If  Wheeler's  orbit  was  retrograde,  or  a  stable  retrograde 
orbit  was  found  nearby,  the  total  AV  required  could  be 
reduced  even  more  as  AV^  for  the  nearby  points  would  be 
substantially  less  due  to  velocity  alignment. 

Since  the  resulting  minimum  AV  transfer  trajectory 
has  a  small  dependence  on  AV^,  the  choice  of  parking 
orbit  is  essentially  an  independent  one.  In  other  words, 
the  parking  orbit  has  very  little  effect  on  what  transfer 
will  have  the  minimum  AV. 

The  amplitude  of  the  total  AV  curve  (Fig  23)  is 
approximately  330  m/sec.  This  is  a  substantial  savings 
and  justifies  the  need  for  a  minimum  AV  transfer.  The 
savings  is  even  more  dramatic  considering  that  most  of  the 
330  m/sec  savings  is  in  AV2  and  fuel  for  AV2  is  payload 
for  AV^. 

The  TOF  curve  is  very  similar  to  the  AV^  curve. 

This  is  due  primarily  to  the  distance  to  the  target  points. 
The  farther  distances  require  a  longer  time  even  though  a 
greater  initial  velocity  was  used  to  depart  the  Earth. 

While  determining  the  minimum  AV  transfer  to  each 
point  two  different  types  of  minimums  were  encountered. 

One  type  was  a  true  minimum,  in  that,  AV  increased  as  you 
moved  away  from  the  minimum  value.  The  transfer  trajectory 
to  point  60  is  an  example  of  this  type.  (See  Table  3.) 
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Total  AV  Versus  Point  Of  Arrival 


Fig  24.  AV.  Versus  Point  Of  Arrival 


Fig  25.  AV~  Versus  Point  Of  Arrival 


Launch  Angle,  <J> ,  Versus  Point  Of  Arrival 


The  other  type  of  minimum  was  one  that  approached  a  limit¬ 
ing  value.  If  a  value  past  the  minimum  was  utilized  it 
would  not  have  sufficient  velocity  to  reach  the  target 
point.  This  type  is  represented  by  the  minimum  transfer 
trajectory  to  point  10.  (See  Table  4) 

Looking  at  the  results  from  the  transfer  to  point  60, 
(Figs  28-31) ,  more  indications  for  the  dominance  of  AVj 
can  be  seen.  Total  AV  (Fig  28)  and  AV (Fig  29)  are 
nearly  identical  in  shape  while  contrasting  sharply  with 
AV^  (Fig  30) .  This  supports  the  conclusion  that  the 
second  burn,  AV^,  is  dominant  in  determining  the  minimum 
AV. 

The  curves  for  point  10  are  dramatically  different  than 
those  of  point  60.  (See  Figs  32-35)  They  do  not  appear  as 
they  would  even  reach  a  minimum,  even  past  the  limiting 
values  for  AV^.  All  of  the  AV  curves  show  an  increasing 
negative  gradient. 

One  of  the  results  of  this  thesis  is  the  algorithm  that 
produces  the  transfer  trajectory.  However,  it  is  not  a 
very  efficient  technique  and  it  does  not  flow  smoothly  or 
easily.  Too  many  educated  adjustments  were  required  because 
the  B  matrix  was  nearly  singular  or  the  required  accuracy 
could  not  be  achieved.  For  these  two  reasons  a  very  good 
initial  guess  is  needed  to  start  the  algorithm.  Without  a 

good  initial  guess  the  path  to  a  solution  is  very  tedious 
and  time  consuming. 
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Fig  28  Total  Delta  V  Versus  Time  of  Flight,  Point  60 


Pig  29  Delta  V1  Versus  Time  of  Flight,  Point  60 
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Table 


Fig  32  Total  Delta  V  Versus  Time  of  Flight,  Point  10 


1 315 


Fig  33  Delta  V,  Versus  Time  of  Flight,  Point  10 


34  Delta  V2  Versus  Time  of  Flight,  Point  10 
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It  is  obvious  from  Fig  23  that  the  minimum  AV  trans¬ 
fer  to  Wheeler's  orbit  does  not  exist  to  point  10.  Assuming 
that  the  total  AV  curve  acts  like  a  quadratic  near  its 
minimum,  a  quadratic  fit  can  be  done  using  points  0,  10 
and  20  to  approximate  the  minimum  value.  An  equation  of  the 
form 

AV  =  A  (pt)2  +  B  (pt)  +  C 

(58) 

would  result.  Solving  the  following  three  equations  will 
find  the  coefficient  for  Eq.  (58) 

4.025413280  -  02A  +  0B  +  C  (Point  0) 

(59) 

3.972123159  =  100A  +  10B  +  C  (Point  10) 

(60) 

3.984135344  -  400A  +  20B  +  C  (Point  20) 

(61) 

Solving  these  and  placing  the  coefficients  in  Eq.  (58) 
results  in 

AV  =  .000326464  (Pt) 2  -  .0086031924  (Pt)  +  4.025413280 

(62) 


To  find  a  minimum 


d(A) 

£TPt) 


.000652928  Pt  -  .0086031924  =  0 
Pt  =  13.156 


75 


(63) 


This  is  approximately  point  13  of  Wheeler's  orbit  or  when 
arrival  time  equals  .816814090  Tu.  Placing  13  in  Eq.  (62) 
gives 

AV  =  3.968828695  Du/Tu  =  3756.98760  m/sec 
This  point  was  not  examined  by  this  thesis. 
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V  Summary  and  Recommendations 


Summary 

A  transfer  trajectory  with  a  minimum  AV  was  found  to 
Wheeler's  orbit.  It  arrived  at  Wheeler's  orbit  at  point  10 
or  when  t  =  .62831853070  Tu.  The  launch  parameters  for 
this  orbit  were  also  established.  From  the  examination  of 
several  orbits  it  was  demonstrated  that  the  minimum  AV 
curve  for  each  orbit  was  very  dependent  upon  AV2.  One  of 
the  major  problems  encountered  was  that  some  of  the  transfer 
trajectories  did  not  have  a  minimum  AV,  but  instead  had  a 
limiting  value  below  which  the  target  point  could  not  be 
reached.  Because  of  this  problem  the  original  algorithm  to 
find  a  minimum  AV  transfer  trajectory  had  to  be  amended. 
The  resulting  algorithm  requires  a  good  initial  guess  for 
it  to  operate  smoothly. 

Recommendations 

The  primary  problem  with  the  algorithm  established  in 
this  thesis  is  its  inability  to  handle  a  singular  B  matrix. 
The  development  of  a  technique  that  deals  with  a  singular 
B  matrix  would  increase  the  time  efficiency  of  calculating 
a  minimum  AV  orbit. 
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The  transfer  trajectory  found  is  not  the  one  with  the 
lowest  AV  possible.  In  the  vicinity  of  target  10  there 
will  be  a  trajectory  with  a  lower  AV.  This  point  was 
approximated  as  point  13.  Further  searching  in  this  area 
would  give  a  more  definite  minimum  AV  trajectory. 

The  results  from  this  thesis  indicate  that  a  retrograde 
orbit  that  is  near  Wheeler's  would  be  very  attractive  in 
minimizing  the  cost  of  a  transfer  to  it.  The  establishment 
of  a  space  colony  or  military  structure  in  a  retrograde 
orbit  would  substantially  lessen  the  transfer  costs  from  the 
Earth.  Further  research  for  the  existance  of  a  retrograde 
orbit  could  prove  to  be  very  beneficial. 

The  techniques  used  in  this  report  should  also  be 
applicable  if  the  originating  position  was  other  than  a  low 
Earth  orbit.  Other  possible  launch  positions  could  be  on 
the  Moon's  surface  or  that  of  a  large  asteroid.  Both  are 
possible  sources  for  raw  materials  for  a  space  colony  or 
military  facility.  Although  the  results  would  be  somewhat 
different,  the  three  launch  locations  could  be  combined  to 
form  a  complete  minimum  transfer  package  to  Wheeler's  orbit 
or  any  other  desired  location. 
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Appendix  A 


Derivation  of  A  Matrix 


The  A  matrix  is  defined  by  Eg.  (36) 


A(t)  * 


•5ji 


1.  4 
1,  4 


(A-l) 


A  4x4  matrix  results  from  (A-l)  with  the  following 
components : 


(A-2) 


(A-4) 


(A-6) 


(A-3) 


(A-5) 


(A-7) 


=  0 


(A-8) 


*  *  2 
(0+a) z 


(1-lOG  - 


=  1.0 


Gu  +  3 (1-u) (X1-M)2G 


(A-9) 
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(A-10) 


+3G(l-y+Xl) 


GMg  SGtXj^+RggCOS  0)  Mg 
+  - - - 


3* 


3  3Gx2(Xl-y) (1-y)  3Gx2y  {l-y+x1) 
VxZ  =  _  -i-  _ 


+  3 GMg (x2“RgEs  in0) (x^+RggCosG) 


(A-ll) 


3x3 

3x3 


*  0 


(A-12) 


3x. 

Fir 


=  2  (0+a ) 


(A-13) 


9*4 

3x3 


3x, 

FxJ 


(A-14) 


«4 

3x, 


G  3G(l-y)(x2)2  3Gyx2 

-  +  - F -  +  - - - 


3 


GMg  3GMg (x2-RggSin0) 
-  +  — 


(0+a)  —  ( 1  — y )  G 

3 


(A-15) 


3*4 


=  -2 (0+a)  (A- 16) 


3*^ 

sir 


=  0 


(A-17) 


4^,  *2»  *3  an<i  *4  are  t^ie  first-order,  non-linear  equations 
of  motion  defined  by  Eqs.  (29)  to  (32) 
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Appendix  B 


Derivation  of  the  B  Matrix 

Utilizing  Fig.  (4),  expressions  for  the  x,  y  coordi¬ 
nates  of  the  satellite  in  its  staging  orbit  can  be  put  in 
terms  of  the  launch  angle,  . 


x  =  r^  cos  <}>  +  y 

(B-l) 

y  =  r^  sin 

(B-2) 

From  these  equations  the  velocity  of  the 

found  by  differentiating  (B-l)  and  (B-2) . 

satellite  is 

•  •  • 

x  =  r^  cos  4>  —  r^  sin  <J>  +  (co  x  r)  term 

(B-3) 

•  •  • 

y  =  r^  sin  <p  +  r^  <j>  cos  <p  +  (a>  x  r)  term 

(B-4) 

The  (w  x  r)  terms  are  the  result  of  the  x_  - 

e  e 

coordinate  rotating  with  respect  to  the  xE  -  y£  system 
which  in  turn  rotates  with  respect  to  the  x-j-  -  yx  system. 


The 

(0) 

x  r) 

term  from  the  x£ 

-  yE 

system  is  on 

the  order 

of 

(rl 

+  Vi) 

and  is  neglected. 

From 

the  inertial 

frame  the 

motion  of  the  satellite  appears  to  be  a  straight  line  and 
the  (cn  x  r)  term  is  neglected  here,  also. 


83 


AD-A079  BB1  AIR  FORCE  INST  OF  TECH  tfRI SHT-PATTERSON  AFB  OH  SCHOO— ETC  F/6  22/3 

DETERMINATION  OF  THE  MINIMJM  DELTA  V  TRANSFER  TRAJECTORY  FROM  A— ETC (II) 
DEC  79  T  E  WEIMER 

UNCLASSIFIED  AFIT/SA-AA/79D-11  ML 


™  1;  ■ 

■ 

■ 

■ 

■ 

■ 

■ 

END 

2  -80 

the  circular 


/ 


*1 


•  • 

Realizing  that  r.  =  0  and  r.4  =  v  , 

•  *  c 

orbit  velocity,  Eq.  (B-3)  and  (B-4)  become 


x  **  -  v  sin  4 
c 

(B-5) 

• 

y  =  vc  cos  <j> 

(B-6) 

Adding  the  impulse  velocity  change  to  these  terms  gives  the 
initial  transfer  velocity  of  the  satellite  at  tg  . 

x  -  -  vc  sin  4  -  Avx  sin  4  =  -  (vc  +  Av j.)  sin  4  (B-7) 

Y  ■  vc  cos  4  +  Avx  cos  <p  -  (vc  +  Av^  cos  4  (B-8) 

Eqs.  (B-l,  2,  7,  8)  provide  the  initial  conditions 
x(tQ) ,  to  begin  the  integration  scheme. 

To  linearize  the  two-point  boundary  value  problem 
Eq.  (38)  is  utilized.  It  relates  errors  at  t^  to  changes 
at  tQ. 

Sxtt^  =  *(t1,  t0)  5x(tQ)  (38) 


•  —i 

- 

5  x 

fix 

5y 

fiy 

• 

s 

* 

fix 

fix 

- 1 

•  >1 
«o 

• 
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However,  the  initial  conditions,  x(tg),  are  related  to 
the  launch  parameters  <j>  and  Av^  by  Eqs.  (B-l,  2,  7,  8)  . 
Utilizing  Eq.  (45) 


5x(t0)  =  C<$L 


(45) 


where 


3x 

3$ 


3x 

3Av^ 


C 


U 

3x 
3  4> 


3x 

3Av^ 


(B-10) 


so. 


— 

•  — 

*  — 

Sx 

5<j> 

6y 

C 

6i 

5Av^ 

(B-ll) 


Combining  (B-9)  and  (B-ll)  produces 


*5x~ 

' 

~6<p  ~ 

6y 

s 

$ 

C 

5* 

5Av1 

t 

(B-12) 
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Since  position  error  is  all  that  we  are  interested  in, 
all  corrections  to  *  and  £  at  t^  are  ignored.  Also, 
calling  the  ($)  (C)  product  by  the  label,  B,  results  in 


r6xi  -  b  i 
Wn  '  KJt0 

[5‘VJ  fco=  E  ]  H 


The  B  matrix  is  defined  as  follows: 


B11  "  *11C11  +  412C21  +  $13C31  +  ®14C41 


B12  *  *11C12  +  $12C22  +  *13C32  +  414C42 


B21  “  *21C11  +  422C21  +  $23C31  +  $24C41 


(B-13) 


(B-14) 


(B-15) 


(B-16) 


B22  *  421C12  +  422C22  *  *23C32  +  424C42 


(B-17) 


where 


'11  3ijr 


-r.^  sin  $ 


(B-18) 


'12  *  5Sv, 


(B-19) 


c2i  ■  n  =  ri cos  4 


(B-20) 


^22-Sfc,-  ° 


(B-21) 
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Appendix  C 


Computer  Program  to  Calculate  Transfer 


Trajectory  to  Wheeler's  Orbit 


This  Appendix  contains  the  computer  routines  used  in 
this  report.  The  computer  language  utilized  was  Fortran 
Extended  Version  IV,  and  all  work  was  accomplished  on  the 
AFIT  CDC  6600  computers.  ODE  was  the  integration  package 
used  in  the  program.  For  the  integration  steps,  one  time 
unit  corresponds  to  4.699939084  days. 
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Theron  E.  Weimer  was  born  8  August  1948  in  Cedar  Rapids, 
Iowa.  He  graduated  from  Thomas  Jefferson  High  School,  Cedar 
Rapids,  Iowa,  June  1966.  The  next  fall  he  entered  the  Uni¬ 
versity  of  Iowa,  Iowa  City,  Iowa.  After  one  year  of  college 
he  received  an  appointment  to  the  United  States  Air  Force 
Academy  and  entered  June  1967.  He  graduated  from  the  USAFA 
June  1971  with  a  Bachelor  of  Science  degree  in  Aeronautical 
Engineering  and  a  regular  commission  as  a  2nd  Lieutenant  in 
the  United  States  Air  Force.  His  initial  assignment  was  to 
Undergraduate  Pilot  Training  (UPT)  at  Williams  Air  Force  Base, 
Arizona.  Upon  graduation  from  UPT  in  June  1972  as  an  Air 
Force  pilot,  he  was  assigned  to  the  316th  Tactical  Airlift 
Wing,  Langley  Air  Force  Base,  Virginia  as  a  C-130E  co-pilot. 
While  there  he  upgraded  to  the  position  of  Aircraft  Commander 
and  also  performed  as  an  airlift  mission  scheduler. 

In  July  1975  he  was  reassigned  to  the  35th  Training 
Squadron,  Reese  Air  Force  Base,  Texas,  as  a  T-37B  instructor 
pilot.  During  the  assignment  at  Reese  he  was  an  assistant 
flight  commander  and  a  Standardization/Evaluation  Flight 
Examiner.  He  also  earned  a  Master  of  Science  degree  in  Bus¬ 
iness  Management  from  the  University  of  Northern  Colorado 
while  at  Reese  Air  Force  Base. 

In  June  1978  he  was  assigned  to  the  Air  Force  Institute 
of  Technology's  resident  School  of  Engineering  at  Wright- 
Patterson  Air  Force  Base,  Ohio  and  began  his  studies  toward 
a  Master  of  Science  degree  in  Astronautical  Engineering. 

Upon  graduation,  he  will  be  assigned  to  the  Combat  Operations 
Center  of  the  Cheyenne  Mountain  Complex,  Colorado  Springs, 
Colorado. 

Permanent  Address:  1615  13th  Street  N.W. 

Cedar  Rapids,  Iowa  52405 
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In  this  report, the  equations  of  motion  for  a  satellite  in 
restricted  four-body  motion  are  derived  and  examined.  An  algorithm 
is  developed  that  produces  a  minimum  delta  v  transfer  trajectory 
from  a  low  Earth  orbit  to  establishment  in  Wheeler's  stable 
periodic  orbit  around  L4.  Both  a  Hohmann  transfer  and  an  infinite 
velocity,  straight-line  transfer  are  examined  as  initial  conditions 
for  integration  of  the  equations  of  motion.  Transfer  trajectories 
are  determined  to  ten  points  around  Wheeler's  orbit  so  that  a. 
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"minimum  delta  v  trajectory  can  be  determined.  Two  types  of 
minimum  delta  v  values  were  found.  One  type  is  a  true  minimum 
value ,  the  other  type  is  a  limiting  value  beyond  which  the 
target  orbit  can  not  be  reached.  Results  sure  presented  on  a 
transfer  trajectory  that  gives  a  minimum  delta  v  from  the  points 
examined.  From  these  results  a  more  accurate  transfer  location 
is  estimated.  _ 
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